Repository navigation
feat(inversion): sparse log det of the split-regularization H on the numpy path (#267 lever 2) - #554
Merged
Merged
Conversation
…numpy path log_det_regularization_matrix_term dense-factorised regularization_matrix_reduced, a (P, P) matrix with O(1) non-zeros per row under the split stencil (~37 ms of the ~300 ms HST Delaunay N=1500 numba CPU call). On the numpy path with the default "cholesky" method it now uses a sparse SuperLU factorisation (scipy splu in symmetric mode, log det = sum(log(diag(U)))), fed by a single mask scan that builds the CSC triple directly: 36.9 ms -> 6.1 ms (6.1x) on the production H (1500x1500, 8.45 non-zeros per row, cond 6.4e12) at 1 thread. It matches the dense Cholesky log det to 7.4e-10 relative on that matrix (-5.7e-6 nats on a log det of 7756.6), which is -2.9e-6 nats and 9.8e-11 relative on the S3 route-b fit's figure_of_merit of 29140.3. At cond 6.4e12 neither fp64 factorisation is the truth to 1e-10 -- numpy's own slogdet differs from the dense Cholesky by up to 2.8e-8 on matrices of this conditioning -- so the difference is the matrix's round-off, not an error introduced on top of it. Two named regime switches keep the dense Cholesky where it is faster, both paid for by the one mask scan: SPARSE_LOG_DET_MIN_PIXELS = 256 (sparse/dense measured at 0.07x for P=9, 0.38x at 128, 1.77x at 256, 10.8x at 1500 -- SuperLU's fixed setup is not worth paying below a few hundred pixels, which also leaves every small fixture, including the numpy/JAX parity fixtures at P=9/16 whose log determinants are pinned across backends at 1e-10, on the historical value), and SPARSE_LOG_DET_MAX_NNZ_PER_ROW = 32 (a kernel scheme's H is coefficient * C^-1 and fully dense, where the sparse route is 7x *slower* at P=1500). The "slogdet" opt-in path and the JAX path are unchanged, as is log_det_curvature_reg_matrix_term; non-PD matrices still raise np.linalg.LinAlgError -- SuperLU's singular RuntimeError is converted and the diagonal pivots are checked, since an LU factorises an indefinite matrix happily where a Cholesky raises -- so the test-mode guard applies unchanged. The stale "uses scipy sparse linear algebra" docstring now describes what runs. Lever 2 of PyAutoLabs/autolens_profiling#267 (stacked on lever 1, PyAutoArray#553). Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Collaborator
Author
|
Paired measurement PR (scripts, JSONs, note section): PyAutoLabs/autolens_profiling#271 — RAL jobs 343353 (CPU A/B, 1.13x, control = lever 1) and 343354 (A100, identity; CPU-only lever). |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Lever 2 of PyAutoLabs/autolens_profiling#267 (phase 3 of the
fixed-lens-light-numba-cpuepic), stacked on lever 1 (#553):AbstractInversion.log_det_regularization_matrix_termno longer dense-factorises the split-regularizationHon the numpy path.After lever 1,
inversion.log_det_regularization_matrix_termwas ~37 ms of a ~300 ms production numba CPU likelihood call (HST Delaunay N=1500,AdaptSplit,InversionImagingSparseNumba, 1 thread): a denseO(P³)Cholesky ofregularization_matrix_reduced, a (1500, 1500) SPD matrix with 8.45 non-zeros per row. The docstring claimed the term "uses scipy sparse linear algebra"; it did not.What changed
inversion_util.log_det_sparse_spd_from(matrix) -> Optional[float]: one pass over the dense matrix builds the CSC triple directly (nocsc_matrix(dense), which costs 13 ms), SuperLU in symmetric mode (permc_spec="MMD_AT_PLUS_A",diag_pivot_thresh=0.0),log det = sum(log(diag(U))). ReturnsNonewhen the dense Cholesky is the faster route.SPARSE_LOG_DET_MAX_NNZ_PER_ROW = 32: geometric-stencilH(neighbour and split schemes, ~5–30 nnz/row, planar-graph structure) goes sparse; kernelH(MaternKerneland siblings, fully dense) stays dense. Real Matern P=1500: dense 34.9 ms vs forced-sparse 335 ms, so the gate matters.SPARSE_LOG_DET_MIN_PIXELS = 256: SuperLU's fixed setup loses below that (sparse/dense 0.38x at P=128, 1.8x at 256, 10.8x at 1500)._log_det_symmetric_from(matrix, sparse=False)offers the matrix to the helper first whensparse=True;log_det_regularization_matrix_termpassessparse=self._xp is np. The"slogdet"opt-in path, the JAX path, andlog_det_curvature_reg_matrix_term(F + λH, dense, lever 3) are byte-identical to lever 1.np.linalg.LinAlgError, exactly asnp.linalg.cholesky, so the test-mode guard applies as before. No fallback, no warning-and-continue.log_det_regularization_matrix_termreplaced with what runs on each backend.API Changes
None in signature or exceptions. Value change on the numpy path only, at round-off level: the sparse LU and the dense Cholesky are two fp64 factorisations of a matrix with condition number 6.4e12 (the
2e-8diagonal jitter against the stencil's smooth null modes), so they agree to 7.4e-10 relative on the log-det term (−5.7e-6 nats) and 9.8e-11 relative on the fit's log evidence (−2.9e-6 nats) at HST Delaunay N=1500 — inside the phase's 1e-9 pin. Below P=256 the value is bit-identical to before (dense path). The existing numpy-vs-JAX parity tests (P=9/16, pinned at 1e-10) are unaffected for that reason; for the record, at that size numpy-splu vs jax-dense differed by 4.3e-10 relative, i.e. a 1e-10 cross-backend pin is only met when both sides run the same algorithm.Measurements
Local (1 thread, real
AdaptSplitH, P=1500): dense Cholesky 36.9 ms → sparse 6.1 ms (of which 2.9 ms is the non-zero scan). RCM + banded Cholesky was tried and rejected (43 ms; the Hilbert ordering leaves a wide band).RAL A/B (control = lever 1
bae9296e, feature = lever 2,gpupartition CPUs-only, 1 thread, phase-2 Leg D cell, route b memo ON):log_det_regularization_matrix_term(ms)bae9296e0b17c2921.13x on the whole call (1.135x on medians); the site falls from 12.4 % to 2.4 % of the call and nothing else moves (the sibling
log_det_curvature_reg_matrix_termis flat at 40.2 → 39.6 ms). The CPU job ran alone on the node (load average at the timed rows 0.49 control, 0.78 feature, no second tenant), so this is a measurement, not a bound. The control reproduces lever 1's feature arm to +1.0 %. Cumulative levers 1+2 against the merge base: 413.3 → 267.4 ms, 1.55x.Witness inside the feature arm on RAL: the sparse route fires exactly once per
figure_of_merit; sparse vs dense-forced log det Δ −9.2e-6 nats (1.2e-9 relative, H eigenvalue ratio 6.5e12);figure_of_meritΔ −4.6e-6 nats (2.3e-10 relative, inside the phase's 1e-9 pin); both gates admit the matrix (P=1500, 8.49 nnz/row); same-matrix timing 6.27 vs 37.72 ms.A100 row (job 343354,
fixed_light.pyDelaunay cell, fp64): identity on the JAX path — all twelvelog_evidence_termsbit-identical between arms. The A100's own dense Cholesky of the same H costs 1.13 ms, 5.6x faster than the CPU sparse route, so this lever is CPU-only by measurement and is deliberately not carried to the GPU.Testing
test_autoarray/inversion/inversion/test_log_det_sparse.py(17 tests): helper vsslogdetat rtol 1e-12 on seeded well-conditioned banded SPD matrices (cond < 1e6 asserted); vs dense Cholesky on realConstantSplit(1e-10) andAdaptSplit(1e-7, with the dense-vs-slogdetspread asserted at the same tolerance so the bound is the matrix's fp64 ambiguity, not a choice) split-stencil H; real kernel H takes the dense route; small H takes the dense route; singular and non-PD matrices raiseLinAlgErroron both routes;"slogdet"path unchanged; the term on a real numpy-path inversion matches its dense value.pytest test_autoarray/inversion -q— 570 passed;pytest test_autoarray/inversion/inversion -q— 229 passed.Tracked by PyAutoLabs/autolens_profiling#267 (lever 2 of 3). Merge order: #553 first.
Heart gate
pyauto-heart readinesswas RED (score 20) at ship time with two release-chain reasons unrelated to this branch — "install verification FAILED (testpypi; checks F)" and "release validation FAILED (stage integrate)" — acknowledged by the human on 2026-09-16 for PR-open only (levers 2 and 3, same reasons as lever 1). Merge stays a human act via/prm.🤖 Generated with Claude Code